System for image reconstruction

ABSTRACT

A method for imaging, including counting quanta of energy emitted into a range of angles from particles of an energy emitter distributed over a volume ( 32 ) in a body ( 30 ), thereby generating a set of counts. The method further includes defining a probability distribution expression that specifies a local concentration of the particles of the energy emitter over the volume as a function of the set of counts, the function being defined in terms of respective coefficients of a plurality of different scales of the local concentration, including at least a first and a second scale. A dependence of the coefficients of the first scale on the coefficients of the second scale is specified, and the local concentration over the volume is computed by applying the probability distribution expression to the set of counts subject to the specified dependence.

CROSS-REFERENCE TO RELATED APPLICATION

This application claims the benefit of U.S. Provisional PatentApplication 60/676,311, filed May 2, 2005, which is incorporated hereinby reference.

FIELD OF THE INVENTION

The present invention relates generally to image processing, andspecifically to imaging using single photon emission sources.

BACKGROUND OF THE INVENTION

Methods for imaging an object emitting photons suffer from a number oflimitations, especially in the case where the object is a human patient,or a section of the patient. The photon emitter is typically aradiopharmaceutical that is injected into, or inhaled by, the patient.In order to minimize deleterious effects on the patient, the quantity ofphoton emitter, the time during which the emitter is active, and thenumber of photons emitted, all need to be minimal. Furthermore,detected, detection systems for the photons typically have minimal or nofocusing ability. Common detection systems may rely on the equivalent ofa fly's eye detection system, wherein the photons are collimated alongtubular channels before detection. As one amongst a number of furthercomplications, the object being imaged is typically three-dimensional,but if a collimator is used the detection system has at best only a veryapproximate way of determining from where, along an axis or within acone of view of the collimator, a detected single photon has beenemitted.

Geometric factors for fly's eye collimators may be calculated. In apaper by Metz et al., entitled “The Geometric transfer functioncomponent for scintillation camera collimators with straight parallelholes,” Phys. Med. Biol., 1980, v. 25, p. 1059-1070, the authors developa method for predicting a geometric transfer function component forconventional scintillation camera collimators. The paper is incorporatedherein by reference. In U.S. Pat. No. 6,943,355, to Shwartz, et al.,whose disclosure is incorporated herein by reference, a method isdescribed that enables a collimator to usefully detect photons fromincident angles exceeding 50.

Thus, geometric factors may be allowed for, given a scan time that issufficiently long. However, the scan time is critical, since the longscan time needed by the imaging system means that the patient is exposedfor the same long time. Unfortunately, raw data produced by operating atshorter times, or with reduced concentrations of radioactivity toproduce the same benefit, has a lower signal to noise ratio than at thelonger times. This acquired raw data is insufficient to enable a goodreconstruction of local concentrations of the radiopharmaceutical, whichis the goal of the imaging system, so that reconstructed images at theshorter times or with reduced concentrations are of extremely lowquality.

Regularization, i.e., reformulation of the raw data, may be used toderive images from the raw data. For example, U.S. Pat. Nos. 5,912,993,6,353,688, and 6,490,374 to Puetter et al., whose disclosures areincorporated herein by reference, describe methods for planar imagereconstruction using “Pixon” elements and bases, which are respectivelydefined as indivisible units of information and sets of possiblefunctions from which the elements are selected.

Wavelet theory may also be used for regularization as described, forexample, by Mallat, in a paper entitled “A theory for multiresolutionsignal decomposition: the wavelet representation,” IEEE Trans. Pat.Anal. Mach. Intell., vol. 11, pp. 674-693, 1989, and by Simoncelli etal., in a paper entitled, “Noise removal via Bayesian wavelet coring,”3rd IEEE Int'l. Conf. Image Processing, Lausanne, Switzerland, September1996, vol. 1, pp. 379-382, IEEE Sig. Proc. Society. Both of these papersassume a Laplacian wavelet distribution. A paper entitled, “Waveletthresholding via a Bayesian approach,” by Abramovich et al., J. R.Statist. Soc. B, vol. 60, pp. 725-749, 1998, describes a variation onthis process, wherein a threshold is set to the wavelet coefficients.The above papers are incorporated herein by reference.

SUMMARY OF THE INVENTION

In embodiments of the present invention, an energy emitter such as aradiochemical is distributed over a specific region, typically withinthe body of a patient. A camera for imaging the region measures numbersof energy quanta from the emitter that are received at the camera. Themeasurements produce a set of counts over multiple angles subtended bythe camera at the specific region. An optimization objective function oflocal concentration of particles of the emitter is defined as a functionof the set of counts and of a number of parameters associated with thecamera, including the geometric configuration of the camera.

The optimization objective function is rewritten as a probabilitydistribution expression. This expression is defined in part in terms ofrespective wavelet coefficients of a plurality of different “scales” ofthe local concentrations at different locations of the specific region.The locations are mapped in multiple dimensions. The term scalecorresponds to a scaling factor applied to location coordinates alongall the dimensions. A dependence between the coefficients of thedifferent scales is defined. A processing unit analyzes the set ofcounts using the expression and the predefined inter-scale dependence,and thus determines values of the local concentrations having thehighest probability in view of the collected counts. The valuescorrespond to a three-dimensional image of local concentration of theenergy emitter within the region. By assuming the predetermineddependence for the scale coefficients, the image produced by theexpression analysis is significantly improved in quality compared toimages from expressions that do not assume a dependence. In cases wherethe energy emitter is a radiochemical, the improvement allows forsignificant reductions in exposure time to the radiochemical and/or theamount of radiochemical injected into the patient.

In one embodiment, the predetermined dependence sets a ratio of thecoefficients of sequential scales to be equal to the expression

$\frac{s^{2}}{( {s + 1} )^{2}},$

wherein s is a cardinality of the wavelet scale.

Typically, the probability distribution expression is defined in termsof a mixture of stationary Gaussian and Poisson distributions relatingthe measured set of counts to an expected set of counts. The expectedset of counts is determined from, inter alia, the geometricconfiguration of the camera. In a disclosed embodiment, for each of theset of counts, the variance of the distribution is equated to themaximum of the count value and a user-determined constant.

In some embodiments of the present invention, the probabilitydistribution expression is defined in terms of the Haar wavelet. Usingthe Haar wavelet enables the wavelet coefficients and their dependenceto be formulated simply, and allows the processing unit to perform therequired analysis in an efficient manner.

In an alternative embodiment, characteristics of elements in thespecific region, such as the motion or viability of the elements, aredetermined by measuring the local concentrations as they change withtime. Values of the measured local concentrations at different times maybe adjusted so that a mean of the measured local concentrationscorresponds to an expected mean. The expected mean is derived from totalcounts over the whole time of measurement.

There is therefore provided, according to an embodiment of the presentinvention, a method for imaging, including:

counting quanta of energy emitted into a range of angles from particlesof an energy emitter distributed over a volume in a body, therebygenerating a set of counts;

defining a probability distribution expression that specifies a localconcentration of the particles of the energy emitter over the volume asa function of the set of counts, the function being defined in terms ofrespective coefficients of a plurality of different scales of the localconcentration, including at least a first and a second scale;

specifying a dependence of the coefficients of the first scale on thecoefficients of the second scale; and

computing the local concentration over the volume by applying theprobability distribution expression to the set of counts subject to thespecified dependence.

Typically, the dependence includes a first cardinality of the firstscale and a second cardinality of the second scale.

The first scale and the second scale may include sequential scales, andthe dependence may be given by:

$\alpha_{s} = \frac{s^{2}}{( {s + n} )^{2}}$

wherein s represents a cardinality of the first scale, 0.5≦n≦1.5, andwherein α_(s) represents the dependence.

In one embodiment an argument of the function comprises a term:

(θ_(i) ^(s)−α_(s)θ_(i) ^(s+1)),

wherein:

i is an index of a location of the particles;

s and s+1 are respective cardinalities of the first and second scales;

θ_(i) ^(s),θ_(i) ^(s+1) represent respective coefficients of the firstand second scales; and

α_(s) is a function of s.

In an alternative embodiment

${f_{s}(z)} = \frac{{{z}_{2}}^{r}}{s^{r}\sigma}$

wherein:

f_(s)(z) represents the function, having an argument z;

s is a cardinality of the first scale;

r is a positive number; and

σ is a constant.

In some embodiments

z=(θ_(i) ^(s)−α_(s)θ_(i) ^(s+1)),

wherein:

i is an index of a location of the particles;

s+1 is a cardinality of the second scale;

θ_(i) ^(s),θ_(i) ^(s+1) represent respective coefficients of the firstand second scales; and

α_(s) is a function of s.

Typically, the probability distribution expression includes a term

−log [P(X)],

wherein

${{- {\log \lbrack {P(X)} \rbrack}} = {\sum\limits_{s}{\sum\limits_{i}{f_{s}( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )}}}},$

and wherein:

X is a vector representing the local concentration;

f_(s)(.) represents the function;

i is an index of a location of the particles;

s and s+1 are respective cardinalities of the first and second scales;

θ_(i) ^(s),θ_(i) ^(s+1) represent respective coefficients of the firstand second scales; and

α_(s) is a function of s.

Alternatively or additionally, the probability distribution expressionincludes a term

−log [P(Y|X)],

wherein

${{- {\log \lbrack {P( {YX} )} \rbrack}} = {\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\Delta_{b}}}}},$

and wherein:

b is an index representing bins receiving the set of counts;

b_(M) is a total number of the bins;

y_(b) is a count received at a bin b;

(HX)_(b) is an expected number of counts for the bin b

Y is a vector representing the set of counts;

X is a vector representing the local concentration; and

Δ_(b) is a variance for the bin b. Typically,

Δ_(b)=max(y _(b) ,B),

wherein B is a constant.

The probability distribution expression may include a term having adistribution chosen from a Poisson distribution, a Laplaciandistribution, and a non-white Gaussian distribution.

Typically, the coefficients include wavelet coefficients of a Haarwavelet.

In a disclosed embodiment, computing the local concentration includesperforming an iteration on the probability distribution expression tofind a maximum value of a probability distribution in the expression.

In some embodiments the local concentration includes a plurality oftime-dependent-local-concentrations including a firsttime-dependent-local-concentration measured at a region of the volumeduring a first time slot and a secondtime-dependent-local-concentration, different from the firsttime-dependent-local-concentration, measured at the region during asecond time slot, and the probability distribution expression includes atime-dependent regularization function relating the firsttime-dependent-local-concentration and the secondtime-dependent-local-concentration.

The time-dependent regularization function typically includes a Fouriertransform of a sequence of time-dependent-local-concentrations at theregion, the sequence including the firsttime-dependent-local-concentration and the secondtime-dependent-local-concentration.

The method typically also includes synchronizing at least one of thefirst time slot and the second time slot in response to a time-dependentsignal received from the body. Typically, computing thetime-dependent-local-concentrations includes iteratively computingpartial sums of the time-dependent-local-concentrations.

In some disclosed embodiments the plurality of different scales includesa number of scales chosen from between three and seven scales, and istypically five scales.

In a further alternative embodiment, the probability distributionexpression includes a conditional expectation.

There is further provided, according to an embodiment of the presentinvention, imaging apparatus, including:

a camera which is arranged to count quanta of energy emitted into arange of angles from particles of an energy emitter distributed over avolume in a body, thereby generating a set of counts; and

a processing unit which is configured to:

define a probability distribution expression that specifies a localconcentration of the particles of the energy emitter over the volume asa function of the set of counts, the function being defined in terms ofrespective coefficients of a plurality of different scales of the localconcentration, including at least a first and a second scale, specify adependence of the coefficients of the first scale on the coefficients ofthe second scale, and

compute the local concentration over the volume by applying theprobability distribution to the set of counts subject to the specifieddependence.

Typically, the local concentration includes a plurality oftime-dependent-local-concentrations including a firsttime-dependent-local-concentration measured at a region of the volumeduring a first time slot and a secondtime-dependent-local-concentration, different from the firsttime-dependent-local-concentration, measured at the region during asecond time slot, and the apparatus includes:

a detector which detects a time-dependent signal from the body and whichconveys the time-dependent signal to the processing unit so as tosynchronize at least one of the first time slot and the second timeslot.

There is further provided, according to an embodiment of the presentinvention, a computer software product for imaging, including acomputer-readable medium having computer program instructions recordedtherein, which instructions, when read by a computer, cause the computerto:

receive counts of quanta of energy emitted into a range of angles fromparticles of an energy emitter distributed over a volume in a body,thereby generating a set of counts,

define a probability distribution expression that specifies a localconcentration of the particles of the energy emitter over the volume asa function of the set of counts, the function being defined in terms ofrespective coefficients of a plurality of different scales of the localconcentration, including at least a first and a second scale,

specify a dependence of the coefficients of the first scale on thecoefficients of the second scale, and

compute the local concentration over the volume by applying theprobability distribution expression to the set of counts subject to thespecified dependence.

The present invention will be more fully understood from the followingdetailed description of the embodiments thereof, taken together with thedrawings, a brief description of which follows.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram of an imaging system, according to anembodiment of the present invention;

FIG. 2 is a schematic diagram illustrating geometric properties of acamera detector in the system of FIG. 1, according to an embodiment ofthe present invention;

FIG. 3 is a flowchart illustrating steps for determining constants inthe system of FIG. 1, according to an embodiment of the presentinvention;

FIG. 4 is a schematic diagram illustrating changes in an image generatedby the system of FIG. 1 with respect to time, according to an embodimentof the present invention;

FIG. 5 is a flowchart illustrating steps for determining gated values ofan image generated in the system of FIG. 1, according to an embodimentof the present invention;

FIGS. 6A and 6B shows schematic images of “cold sphere” phantoms; and

FIGS. 7A and 7B shows schematic images of “hot sphere” phantoms.

DETAILED DESCRIPTION OF EMBODIMENTS

Embodiments of the present invention provide a method of imagereconstruction for images that are generated by particles of energyemitters and detected by detectors at different locations. For example,in different aspects of the present invention, the energy emitterparticles may comprise a radiochemical which emits electromagnetic (EM)energy by radioactive decay; an excited atom, such as is formed in amagnetic resonance environment, which emits EM energy by changing itsenergy state; or a molecule which emits sound energy in response toultrasonic waves impinging on the molecule. The emitted energy may inall cases be assumed to be in the form of photons or phonons. In thespecification and in the claims, the terms quanta or energy quanta areassumed to encompass both photons and phonons. In the specification andin the claims, the term energy emitter is assumed to encompass bothactive emitters such as the radiochemical exemplified above, and passiveemitters, such as the molecule emitting sound energy as reflectedultrasound. The terms “camera” or “camera detector” refer to any sort ofdetection system suitable for measuring the emission of the quanta inquestion.

It will be understood that while the description below refers in generalto a radiochemical as the energy emitter, the principles of the presentinvention are applicable, mutatis mutandis, to energy emitters of othertypes.

The problem of reconstructing local concentrations of particles of aradiochemical from the raw data generated by a camera is an example ofan ill-posed inverse problem, in the sense that there may not be asingle solution, and even if a single solution exists, it may beunstable with respect to small perturbations of the raw data results.

Methods for improving the quality of the results generated by thephysical detection system described above have been sought. Some ofthese methods are based on the application of Bayes' Theorem:

$\begin{matrix}{{P( I \middle| D )} = \frac{{P( D \middle| I )} \cdot {P(I)}}{P(D)}} & (1)\end{matrix}$

where D represents the raw data, I represents the image, P(I|D) is theprobability of I given D, and P(I) is the probability of I.

An expression (2) derived from equation (1) is:

P(I|D)∝P(D|I)·P(I)  (2)

where the first term on the right side of the proportionality, measuringthe likelihood of D given I, may be referred to as the likelihood term,and the second term may be referred to as the image prior, expressing ana priori expectation of a particular image.

The Bayesian approach to solving ill-posed inverse problems was firstdescribed by Jaynes in his paper “Prior Information and Ambiguity inInverse Problems,” Inverse Problems, D. W. McLaughlin (ed.), Am. Math.Soc., SIAM-AMS Proceedings vol. 14, 1984. Bayesian image reconstructionwas first used by Geman et al. in their paper “Stochastic relaxation,Gibbs distributions, and the Bayesian restoration of images,” IEEETrans. PAMI, PAMI-6(6), November 1984. Both papers are incorporatedherein by reference.

Earlier approaches at regularization of raw data use wavelets, andtypically assume that results follow a Gaussian distribution. However,more recent attempts at regularization using wavelet theory have assumedother distributions.

The simplest wavelet, the Haar wavelet, is defined as:

$\begin{matrix}{{h(x)} = \{ \begin{matrix}1 & {0 \leq {\times {< \frac{1}{2}}}} \\{- 1} & {\frac{1}{2} \leq {\times {< 1}}} \\0 & {otherwise}\end{matrix} } & (3)\end{matrix}$

While being the simplest wavelet, the Haar wavelet has the disadvantagethat it is discontinuous, and therefore not differentiable at thediscontinuities.

Reference is now made to FIG. 1, which is a schematic diagram of animaging system 20, according to an embodiment of the present invention.System 20 comprises a camera 22 coupled to a processing unit (PU) 24,which processes raw data received from the camera. PU 24 includes aprocessor 26 and a memory 28, wherein are stored software instructionsfor storing and processing the raw data. Processor 26 may be ageneral-purpose computer or may comprise hard-wired or programmablelogic. In the event that processor 26 comprises a general-purposecomputer the software may be provided as a computer software product ina tangible form on a computer-readable medium such as a CD-ROM, or as anelectronic data transmission, or as a mixture of both forms.

In operating system 20, an energy emitter, herein assumed to be aradiochemical formed from a radioisotope such as ^(99m)Tc or ¹²³I, isinjected into a patient 30, and gamma ray quanta emitted by the decay ofthe radioisotope are detected by camera 22. The raw data recorded bycamera 22 is reconstructed by PU 24 to determine a three-dimensional mapof the concentrations of the radioisotope, and hence of theradiochemical, in a section 32 of the body of patient 30 that is scannedby camera 22. PU 24 may display the three-dimensional map as atwo-dimensional image on a graphic display unit 23, for viewing by anoperator 21 of system 20.

Camera 22 comprises a detector 34 which is divided into atwo-dimensional array 40 of pixels. Detector 34 typically comprises oneor more scintillators followed by respective photo-multipliers.Alternative detectors which do not use photo-multipliers, such assolid-state detectors, are known in the art, and it will be appreciatedthat such detectors are within the scope of the present invention. Inthe following explanation, detector 34, by way of example, is assumed touse scintillators. Optionally, each scintillator may be configured tocount gamma ray photons only within a specific energy range. The rangemay be defined so as to exclude energies of photons that may haveundergone scattering. (Scattered photons typically have energies lowerthan unscattered photons.)

A collimating grid 36, before detector 34 comprises a plurality ofcollimators 38, typically in the form of closely packed hexagons,although the collimators may be any other type known in the art, such ascylinders or cones. Typically, array 40 of pixels is configuredindependently of collimators 38. Optionally, processing unit 24 usescollimators 38 to divide detector 34 into the array of detector pixels,and may allocate a group of contiguous collimators to one pixel. By wayof example, array 40 is assumed to comprise a square array of 64×64pixels, but it will be understood that the number of pixels, and theirgeometric arrangement and shape, may be different from that exemplifiedhere, and may comprise any convenient number and geometry.

There are a number of methods, known in the art, by which camera 22acquires its raw data. For example, camera 22 may be held stationary inrelation to patient 30, may be continuously moved in relation to thepatient, or the camera may be moved from one position to another, beingheld stationary at each position. Herein, by way of example, an operatorof system 20 is assumed to rotate camera 22 in a 180° arc of a circle,perpendicular to an axis 42 of patient 30, and the camera is heldstationary at sixty positions separated by 3° on the arc. The camerarecords its raw data at the sixty different positions. Camera 22 thusgenerates 64×64×60=245,760 values of raw data. In an alternativeembodiment, the camera records 64 sets of data, generating 64³=262,144values of raw data. However, it will be understood that the total numberof values of raw data depends on the number of pixels defined for camera22, as well as on the number of positions of the camera, and thatsubstantially all such numbers are within the scope of the presentinvention. It will also be understood that the raw data may be acquiredby substantially any method known in the art, and that all such methodsare assumed to be comprised within the scope of the present invention.

PU 24 receives the raw data, and processes it to find concentrations ofthe radiochemical in section 32. The raw data is herein represented as avector Y={y_(b)}, where b is a raw data bin index from 1 to b_(M), andwhere b_(M) represents the total number of values of raw data.Optionally, raw data Y may be pre-processed by PU 24 before beinganalyzed as described below. Such pre-processing is well known in theart, and typically includes rearrangement of actual numbers of countsdetected by detector 34 into bins b, noise reduction, and/or detectordistortion compensation, etc.

In the computation that follows, section 32 is taken to comprise a cubeformed of 64³, 262,144, voxels. Each voxel is a cube having apre-determined edge length, which is assumed by way of example to beapproximately 0.68 cm. It will be appreciated, however, that section 32and/or the voxels forming the section may have a geometric shapedifferent from a cube; for example the section and/or the voxels may bebox-shaped, with the edges of the box having different values; it willalso be appreciated that the number of voxels within section 32, and thesize of the voxels, may be different from the values exemplified herein.

The concentration of particles of radiochemical in section 32 is hereinrepresented as a vector X={x_(i)}, where i is an index from 1 to v_(T),and where v_(T) represents the total number of voxels. It will beunderstood that elements x_(i) of vector X are distributed in threespatial dimensions and that i determines a location of x_(i).

Optionally, a detector 44 is coupled to patient 30, typically for thepurpose of allowing PU 24 to correlate patient parameters measured bythe detector with raw data Y. An example of the use of such a detectoris described below with reference to FIG. 4.

FIG. 2 is a schematic diagram illustrating geometric properties ofcamera 22, according to an embodiment of the present invention. Asection 50 of array 40, when camera 22 is in a given position,illustrates four bins b−1, b, b+1, b+2. Each bin corresponds to one ofthe pixels of array 40, for the given position of the camera. A portion60 of section 32 includes four voxels i−1, i, i+1, i+2 of the section.Each voxel radiates gamma ray photons in proportion to the concentrationof radiochemical in the voxel, but the number of photons detected byeach of the bins, y_(b−1), y_(b), y_(b+1), y_(b+2), depends on factorsother than the concentrations, including geometric factors. For example,solid angles A_(i−1,b), A_(i,b), A_(i+1,b), A_(i+2,b) are subtended byeach of the voxels i−1, i, i+1, i+2 at bin b, and the number of photonsdetected by bin b is a function of the solid angles.

As will be apparent to those skilled in the art, other factors which maydetermine the number of photons detected by a specific bin include, butare not limited to, the following:

-   -   An effective area of the bin, which typically takes into account        a blocking effect of collimators. Geometric effects of the        collimators may be allowed for, as described in the Background        of the Invention.    -   A probability that a photon may penetrate septa comprised in        collimators of the bin.    -   Attenuation and/or scattering of photons due to interaction of        the photons with media between the bin and a radiating voxel.    -   Energy of the emitted photons, since the energy may affect the        values of effective area, probability, attenuation and        scattering listed above, as well as a response function of        detector 34.

A vector Y={ y_(b) } represents an expected number of photons detected,in the absence of noise, in each bin b due to the factors describedabove. In embodiments of the present invention, we assume:

Y=HX  (4)

where H is a matrix having coefficients h_(b,i).

Coefficients h_(b,i) take into account part or all of the factorsdescribed above. Following from equation (4) the notation (HX)_(b) isalso used for y_(b) .

For specific data Y, we assume that the required value of each x_(i) isgiven by:

{x _(i)}_(i=1) ^(v) ^(T) =argmax[P(X|Y)]  (5)

where P(X|Y) is the probability of the occurrence of X given Y.

In other words, the values of x_(i) correspond to the vector X that hasa maximum probability of occurring given the observed data Y.

From equation (1), P(X|Y) may be rewritten:

$\begin{matrix}{{P( X \middle| Y )} = \frac{{P( Y \middle| X )} \cdot {P(X)}}{P(Y)}} & (6)\end{matrix}$

where P(X) and P(Y) are the respective probabilities of X and Yoccurring.

By taking logarithms, we derive the following expressions from equation(6):

log [P(X|Y)]=log [P(Y|X)]+log [P(X)]−log [P(Y)]  (7)

Since Y is the data received by PU 24, log [P(Y)] is a constant, which,as is explained below, may be assumed to be 0 without loss ofgenerality. Applying this to equation (7), and multiplying by −1, gives:

log [P(X|Y)]=−log [P(Y|X)]−log [P(X)]  (8)

Referring back to equation (5), finding the vector X having a maximumprobability of occurring corresponds to finding the vector X that causes−log [P(X|Y)] to be a minimum. The terms on the right side of equation(8) are termed the likelihood term, −log [P(Y|X)], and the penalty term,−log [P(X)], and are analyzed separately below.

Likelihood Term, −log [P(Y|X)]

In an embodiment of the present invention, the values of y_(b) areassumed to form an approximately Gaussian distribution, so that:

$\begin{matrix}{{P( y_{b} \middle| X )} = ^{\frac{- {({y_{b} - {({HX})}_{b}})}^{2}}{2\Delta_{b}}}} & (9)\end{matrix}$

where Δ_(b) is the variance of y_(b).

From equation (9), assuming the terms are independent, the expressionfor the likelihood term becomes:

$\begin{matrix}{{- {\log \lbrack {P( Y \middle| X )} \rbrack}} = {\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\Delta_{b}}}}} & (10)\end{matrix}$

In the present embodiment,

Δ_(b)=max(y _(b) ,B)  (11)

where B is a predetermined value.

Typically, the value of B is determined according to the actual valuesof y_(b), and is a predetermined constant of the order of 200. The valueof B is set to ensure that low values of y_(b) do not unduly skew thevalue of the likelihood term. A method for determining B is described inreference to FIG. 3 below.

From equation (11), the expression for the likelihood term for thedisclosed embodiment may be rewritten as:

$\begin{matrix}{{- {\log \lbrack {P( Y \middle| X )} \rbrack}} = {\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\max \; ( {y_{b},B} )}}}} & (12)\end{matrix}$

It will be appreciated that expression (12) is one specific example ofthe derivation of the likelihood term, and that the term may be derivedby a number of related methods, giving correspondingly differentexpressions. Examples of some alternative methods for deriving thelikelihood term are as follows:

-   -   The variance Δ_(b) could be defined differently from        equation (8) while still including a lower bound B and having a        value that varies according to the bins. Examples of different        expressions for Δ_(b) are √{square root over (y_(b) ²+B²)} or        (y_(b)+B); other such expressions will be apparent to those        skilled in the art.    -   Rather than having a different variance for each bin b, the        variance may be a value Δ that is the same for all bins.    -   Equation (9) assumes an approximately Gaussian distribution;        other distributions which may be used include a Poisson        distribution, a mixture of a Gaussian and a Poisson        distribution, a Laplacian distribution, or a non-white Gaussian        distribution. In the latter case, the expected value (HX)_(b)        may be assumed to be weighted according to expected values from        neighboring bins, such as (HX)_(b−1) and/or (HX)_(b+1).

Other derivations for the likelihood term will be apparent to thosehaving skill in the art, and all such derivations are assumed to becomprised within the scope of the present invention.

Penalty Term, −log [P(X)]

The inventors have found that writing the penalty term in terms ofwavelet coefficients, where the coefficients at different scales are setto be dependent on each other, results in reconstruction of generallyclearer images than are produced by prior art systems. A specificexample of such a dependency is described in more detail hereinbelow.

A general expression for the penalty term is defined as:

$\begin{matrix}{{- {\log \lbrack {P(X)} \rbrack}} = {\sum\limits_{s}{\sum\limits_{i \in Z^{3}}{f_{s}( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )}}}} & (13)\end{matrix}$

where f_(s)(.) is a function;

s are positive integers representing a cardinality of the scale of thewavelet;

i correspond to the indices, defined above, for vector X; herein i≡(n₁,n₂, n₃)εZ³ are ordered triples of integers n₁, n₂, n₃, identifying thevoxels of X;

θ_(i) ^(s),θ_(i) ^(s+1) are wavelet coefficients; and

α_(s), where α_(s)εR|α_(s)>0, is a term relating the dependency ofsequential wavelet coefficients.

Expression (13) is evaluated in three dimensions (the three spatialdimensions of vector X), so that the ordered triples i, which havevalues from (1, 1, 1) to (64, 64, 64), cover these spatial dimensions.Thus, θ_(i) ^(s) is a three dimensional vector which may be written:

θ_(i) ^(s)={θ_(i,d) ^(s),d=1,2,3}  (14)

The wavelet used herein, by way of example, is the Haar wavelet, definedin equation (3). The properties of the Haar wavelet allow us to rewritethe elements of wavelet coefficients θ_(i) ^(s) in terms of elementsx_(i) as follows:

$\begin{matrix}{\theta_{i,d}^{s} = {{\sum\limits_{m = 0}^{s - 1}x_{i + {m \cdot j_{d}}}} - {\sum\limits_{m = {- s}}^{- 1}x_{i + {m \cdot j_{d}}}}}} & (15)\end{matrix}$

where {j_(d), d=1, 2, 3} is an orthonormal standard basis in Z³.

In this embodiment, the following limitations are placed on theexpression for the penalty term:

1. The dependency term α_(s) is defined in terms of the cardinality ofthe scale:

$\begin{matrix}{\alpha_{s} = \frac{s^{2}}{( {s + n} )^{2}}} & (16)\end{matrix}$

where 0.5≦n≦1.5. Herein it is assumed that n=1.

The argument of function f_(s)(.) (equation 13) is (θ_(i)^(s)−α_(s)θ_(i) ^(s+1)), and expressions for θ_(i,d) ^(s) and α_(s), asdefined by equations (15) and (16), enable us to evaluate the argumentfor predefined functions of x.

2. The function f_(s)(.) is assumed to be of the form:

$\begin{matrix}{{f_{s}(z)} = \frac{{{z}_{2}}^{r}}{s^{r}\sigma}} & (17)\end{matrix}$

Thus,

$\begin{matrix}{{f_{s}( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )} = \frac{{{{\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}}}_{2}}^{r}}{s^{r}\sigma}} & (18)\end{matrix}$

where σ is a constant, and r>0. Herein, r is assumed equal to 1.

The Penalty term may thus be written:

$\begin{matrix}{{- {\log \lbrack {P(X)} \rbrack}} = {\frac{1}{\sigma}{\sum\limits_{s}{\sum\limits_{i \in Z^{3}}\frac{{( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )}_{2}}{s}}}}} & (19)\end{matrix}$

Returning to equation (8), from equations (12) and (19) the expressionfor −log [P(X|Y)] may be rewritten:

$\begin{matrix}{{- {\log \lbrack {P( X \middle| Y )} \rbrack}} = {{\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\max \; ( {y_{b},B} )}}} + {\frac{1}{\sigma}{\sum\limits_{s}{\sum\limits_{i \in Z^{3}}\frac{{( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )}_{2}}{s}}}}}} & (20)\end{matrix}$

Inspection of equation (20) shows that the value of −log [P(X|Y)]depends on two factors, B and σ, which may be set by experimentation. Anexemplary method for determining B and σ is described below withreference to FIG. 3. From equation (20), the value of −log [P(X|Y)] alsodepends, inter alia, on the number of scales s summed in the Penaltyterm. The inventors have found that summing values of s between 1 and 5,as shown in equation (21) below, gives a satisfactory result. Otherranges of values for s include values between 1 and 3, and valuesbetween 1 and 7. However, those skilled in the art will be able toevaluate alternative ranges for s without undue experimentation, and allsuch ranges are assumed to be within the scope of the present invention.

If scales s are summed between 1 and 5, equation (20) may be rewrittenas the expression:

$\begin{matrix}{{- {\log \lbrack {P( X \middle| Y )} \rbrack}} = {{\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\max \; ( {y_{b},B} )}}} + {\frac{1}{\sigma}{\underset{s = 1}{\sum\limits^{s = 5}}{\sum\limits_{i \in Z^{3}}\frac{{( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )}_{2}}{s}}}}}} & (21)\end{matrix}$

Embodiments of the present invention determine values of x_(i) thatminimize −log [P(X|Y)] and thus maximize P(X|Y). (Returning briefly toequation (7), since the term log [P(Y)] of equation (7) is a constant,and since we find a minimum of −log [P(X|Y)], the term log [P(Y)] may beassumed to be 0 without loss of generality.)

The determination of the minimum may be performed using any suitablealgorithm known in the art for finding minima of functions of the typeexemplified by equation (20), including both iterative and non-iterativealgorithms.

One such iterative algorithm, which the inventors have found performsefficiently, is a limited-memory quasi-Newton algorithm. The algorithmis described in “A limited memory algorithm for bound constrainedoptimization,” SIAM Journal of Scientific and Statistical Computing, 16,5, pp. 1190-1208, 1995; and “Representations of quasi-Newton matricesand their use in limited memory methods,” Mathematical Programming, 63,4, pp. 129-156, 1994, both papers by Byrd et al., whose disclosures areincorporated herein by reference. A part of the algorithm is alsodescribed in “Line search algorithms with guaranteed sufficientdecrease,” ACM Transactions on Mathematical Software 20 (1994), no. 3,pp. 286-307, by Mor'e et al., whose disclosure is incorporated herein byreference. In one embodiment of the present invention, the algorithm isrun a preset number of times before stopping. Other suitable stoppingcriteria, such as evaluating a rate of convergence, will be apparent tothose of skill in the art.

The algorithm that determines the minimum of −log [P(X|Y)] returns avalue for X, corresponding to a three-dimensional image of section 32.Hereinbelow, the algorithm used is termed the minimum-finding algorithm.Typically, after using the minimum-finding algorithm, system 20 presentsthe three-dimensional image as two-dimensional slices on unit 23.Optionally, prior to presentation on unit 23, processing unit 24 maypost-process the values of X derived from the minimum-finding algorithm.

In an alternative embodiment of the present invention, a maximum for aconditional expectation E(X|Y), is found. E(X|Y) is defined as follows:

$\begin{matrix}{{E( X \middle| Y )} = \frac{\int_{{x_{j} \geq 0},{j \in Z^{3}}}{{X \cdot {P( X \middle| Y )}}\ {X}}}{\int_{{x_{j} \geq 0},{j \in Z^{3}}}{{P( X \middle| Y )}\ {X}}}} & (22)\end{matrix}$

where P(X|Y) may be found from the expression for −log [P(X|Y)] given inequation (21).

Methods used to find the maximum of E(X|Y) are generally similar tothose for finding the maximum for P(X|Y), described herein.

FIG. 3 is a flowchart 150 illustrating steps for determining σ and B,used in equation (21), according to an embodiment of the presentinvention. Before implementing the minimum-finding algorithm, theconstants σ and B are determined, typically by a method similar toflowchart 150. In an initial step 152 of the flowchart, a database ofraw data is assembled. Typically the database comprises {y_(b)} valuestaken from patients. In one implementation of flowchart 150, theinventors used values of {y_(b)} from 26 patients. Alternatively oradditionally, the database comprises {y_(b)} values from a phantom,which correspond in quality, e.g., count number, to patient data.

In a step 154, a first image is determined from {y_(b)} by applying theminimum-finding algorithm to equation (21), using initial values for σand B.

In a comparison step 156, the first image is compared with an expectedimage. For data such as that from a phantom, the expected image may bededuced from the known geometric distribution of the radiochemical;alternatively, the expected image may be determined by using a largequantity of radiochemical and/or a long exposure time to improve thequality of the results presented to PU 24. For data such as those frompatients, one or more expert physicians may generate the expected image.If the comparison is satisfactory, flowchart 150 ends. If the comparisonis unsatisfactory, the flowchart returns to step 154, and σ and/or B arechanged.

Steps 154 and 156 reiterate until a satisfactory comparison is achievedin step 156.

The inventors have found that in one embodiment of the presentinvention, flowchart 150 generated values for σ and B equal to 1000 and200 respectively.

FIG. 4 is a schematic diagram illustrating changes in X with respect totime t, according to an embodiment of the present invention. Thedescription above does not explicitly incorporate time in finding valuesof X, although in fact X typically does vary with time due to the decayof the radiochemical itself. In addition, depending on where section 32is located in the body of patient 30, and/or of the malady being imaged,there may be a variation of X with time due to, inter alia, movement ofparts of the body. This variation may be periodic, non-periodic, or acombination of the two. For example, if section 32 comprises the heart,there is a periodic change in X due to the heart beating. If section 32comprises a blood vessel, there is movement of blood through the vessel,and the movement may comprise a combination of periodic motion andtranslation of the blood. As a further example, over time, necrotictissue may display differently compared to viable tissue. FIG. 4illustrates the changes of X with time by assuming that there are anumber of discrete values of X, represented by X(t), at time slots t=0,1, 2, . . . . Each X(0), X(1), X(2), . . . , generates respectiveresults {y_(b)}₀, {y_(b)}₁, {y_(b)}₂, . . . . Herein X(0), X(1), X(2), .. . , are also termed gated values of X, and {y_(b)}₀, {y_(b)}₁,{y_(b)}₂, . . . are also termed gated values of Y.

In one embodiment of the present invention, the values of X(0), X(1),X(2), . . . are determined substantially independently of each other, sothat essentially any relation between the different values is ignored.The determinations are by applying the minimum-finding algorithm to eachset of values {y_(b)}₀, {y_(b)}₁, {y_(b)}₂, . . . .

In an alternative embodiment of the present invention, sequential valuesX(0), X(1), X(2), . . . are assumed to have a relation to each other.The following description uses the periodic motion of a beating heart toexemplify the change of X with time t. Those skilled in the art will beable to adapt the description, mutatis mutandis, for other periodic andnon-periodic changes of X with time.

Returning to FIG. 1, detector 44 comprises an electrocardiograph (ECG)sensor, and PU 24 uses an output of the ECG to synchronize counts{y_(b)} with different stages of the beating heart so as to cover acomplete heat beat, i.e., between end-systolic and end-diastolic states.PU 24 uses the ECG output to divide counts {y_(b)} into N time slots{t=0, 1, 2, . . . , N−1}, where N is a positive integer, and PU 24 sumsthe values of the counts for each specific bin b to give N differentvalues {y_(b0), y_(b1), . . . y_(b(N−1))}. In an embodiment of thepresent invention, N may be chosen to be equal to 2^(n), where n is apositive integer, although any other convenient value of N, for example24, may be used.

Taking account of the time dependence of X, equation (20) may berewritten:

$\begin{matrix}{{- {\log \lbrack {P( {XY} )} \rbrack}} = {{\sum\limits_{t = 0}^{N - 1}\begin{pmatrix}{{\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\max ( {y_{b},B} )}}} +} \\{\frac{1}{\sigma}{\sum\limits_{s}{\sum\limits_{i}\frac{{( {\theta_{i}^{2} - {\alpha_{s}\theta_{i}^{s + 1}}} )}_{2}}{s}}}}\end{pmatrix}} + {\sum\limits_{i}{W( \{ {{x_{i\;}(t)},{t = 0},1,{{\ldots \mspace{14mu} N} - 1}} \} )}}}} & (23)\end{matrix}$

where the expressions

$\frac{{( {\theta_{i}^{2} - {\alpha_{s}\theta_{i}^{s + 1}}} )}_{2}}{s}$and$\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\max ( {y_{b} - B} )}$

are evaluated separately for each time slot t;

{x_(i)(t), t=0, 1, 2, . . . N−1} are the sets of x_(i) at time slots t;and

W(.) is a suitable one-dimensional regularization function.

In one embodiment of the present invention, W₍ ₎ is expressed in termsof a Fourier transform of x_(i)(t), and is defined as:

$\begin{matrix}{{W( \{ {{x_{i}(t)},{t = 0},1,\ldots \mspace{14mu},{N - 1}} \} )} = {\sum\limits_{k = 0}^{N - 1}{{\overset{\sim}{g_{i}}( \omega_{k} )}{{\overset{\sim}{x_{i}}( \omega_{k} )}}^{2}}}} & (24)\end{matrix}$

where {{tilde over (x)}_(i)(ω_(k)), k=0, 1, . . . , N−1} is the1-dimensional Fourier transform of the sequence {x_(i)(t), t=0, 1, . . ., N−1};

${\omega_{k} = \frac{2\pi \; k}{N}};{and}$

{tilde over (g)}_(i)(ω_(k)) are N filtering constants (for every i),typically determined as described above with respect to FIG. 3 for B, σ.

In an alternative embodiment, W(.) is a one-dimensional expression basedon equation (19):

$\begin{matrix}{{W( \{ {{x_{i}(t)},{t = 0},1,{{\ldots \mspace{14mu} N} - 1}} \} )} = {\frac{1}{\sigma_{t}}{\sum\limits_{s = 1}^{N/2}{\sum\limits_{k = 0}^{N - 1}\frac{{\begin{pmatrix}{{\theta_{i}^{s}(k)} -} \\{\alpha_{s}{\theta_{i}^{s + 1}(k)}}\end{pmatrix}}_{2}}{s}}}}} & (25)\end{matrix}$

where N is assumed to be even;

${\alpha_{s} = \frac{s^{2}}{( {s + 1} )^{2}}};$ and${\theta_{i}^{2}(k)} = {{\sum\limits_{m = 0}^{s - 1}{x_{i}( {k + m} )}} - {\sum\limits_{m = {- s}}^{- 1}{{x_{i}( {k + m} )}.}}}$

Further alternatively, W(.) may be defined based on the assumption thatresults follow a Laplacian distribution similar to that described in thepapers of Mallat and Simoncelli et al. referenced in the Background ofthe Invention. For example,

$\begin{matrix}{{{W( \{ {{x_{i\;}(t)},{t = 0},1,\ldots \mspace{14mu},{N - 1}} \} )} = {\frac{1}{\sigma_{t}}{\sum\limits_{k = 0}^{N - 1}{{{x_{i}(k)} - {x_{i\;}( {k + 1} )}}}^{r_{t}}}}}{or}{{W( \{ {{x_{i}(t)},{t = 0},1,\ldots \mspace{14mu},{N - 1}} \} )} = {\frac{1}{\sigma_{t}}{\sum\limits_{k = 0}^{N - 1}{{{x_{i\;}( {k - 1} )} - {2{x_{i\;}(k)}} + {x_{i\;}( {k + 1} )}}}^{r_{t}}}}}} & (26)\end{matrix}$

In equations (25) and (26), since {x_(i)} are assumed to be periodic andare measured in N repeating time slots, x_(i)(k) is to be considered asx_(i)(k(mod N)). Also in the equations, σ_(t) and r_(t) are constantswhich are typically determined as described above with respect to FIG.3.

Other expressions for W(.) will be apparent to those skilled in the art,and are comprised within the scope of the present invention.

FIG. 5 is a flowchart 200 showing steps for finding gated values of X,assuming W(.) is of the form given by equation (24), according to anembodiment of the present invention.

From equation (21), a function Φ^(3D)(X,Y,B,σ) is defined:

$\begin{matrix}{{\Phi^{3D}( {X,Y,B,\sigma} )} \equiv ( {{- {\log \lbrack {P( {XY} )} \rbrack}} = {{\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\max ( {y_{b},B} )}}} + {\frac{1}{\sigma}{\sum\limits_{s = 1}^{s = 5}{\sum\limits_{i \in Z^{3}}{\frac{{( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )}_{2}}{s}.}}}}}} } & (27)\end{matrix}$

In a first step 202, the following partial sums of the gated values ofX, and partial sums of gated data Y, are defined:

$\begin{matrix}{{{X^{(1)}(\tau)} \equiv {\sum\limits_{\tau^{\prime} = 0}^{2^{n - 1} - 1}{X( {\tau + {2^{1}\tau^{\prime}}} )}}}{{Y^{(1)}(\tau)} \equiv {\sum\limits_{\tau^{\prime} = 0}^{2^{n - 1_{- 1}}}{Y( {\tau + {2^{1}\tau^{\prime}}} )}}}} & (28)\end{matrix}$

where 0≦1≦n, 0≦τ<2¹, n, 1, and τ are 0 or positive integers. Herein, byway of example, n is assumed equal to 3.

The parameter 1 in equations (28) corresponds to a level of partial sum.Thus, for n=3, the highest value of 1 is 3 and each partial sum is oneterm corresponding to a gated value, e.g., X⁽³⁾(2)≡X(2). For 1=1, thereare four terms in each partial sum of gated values, e.g.,X⁽¹⁾(0)≡X(0)+X(2)+X(4)+X(6).

In a second step 204, a “dc” value for X, X_(dc), corresponding to thevalue of X calculated from the total values in each bin b of y_(b), iscomputed. In addition, an optimal value {circumflex over (X)}⁽⁰⁾(0) isdefined.

$\begin{matrix}{{X_{d\; c} = {\underset{\{{{0 \leq x_{i}^{\prime}},{i \in Z^{D}}}\}}{\arg \mspace{11mu} \min}{\Phi^{3D}( {{X^{\prime};{Y^{(0)}(0)};B},\sigma} )}}}{{{\hat{X}}^{(0)}(0)} \equiv X_{d\; c}}} & (29)\end{matrix}$

It will be understood that step 204 corresponds to evaluating equation(21).

In a third step 206, optimal values of {circumflex over (X)}⁽¹⁾(τ),0≦τ<2¹⁻¹, are evaluated recursively, for 1=1, 2, 3, as follows:

$\begin{matrix}{{\begin{Bmatrix}{{{\hat{X}}^{(1)}(\tau)},} \\{0 \leq \tau < 2^{1 - 1}}\end{Bmatrix} = {\underset{\{\begin{matrix}{{0 \leq {x_{i}^{(1)}{(\tau)}}\; \leq {{\hat{x}}_{i}^{({1 - 1})}{(\tau)}}},{i \in Z^{3}},} \\{0 \leq \tau < 2^{1 - 1}}\end{matrix}\}}{\arg \mspace{11mu} \min}{\Phi^{(1)}( {\begin{Bmatrix}{{X^{(1)}(\tau)},} \\{0 \leq \tau < 2^{1 - 1}}\end{Bmatrix};\begin{matrix}{Y^{(1)},} \\{\hat{X}}^{({1 - 1})}\end{matrix}} )}}}{where}} & (30) \\{{{{\Phi^{(1)}( {\begin{Bmatrix}{{X^{(1)}(\tau)},} \\{0 \leq \tau < 2^{1 - 1}}\end{Bmatrix};\begin{matrix}{Y^{(1)},} \\{\hat{X}}^{({1 - 1})}\end{matrix}} )} \equiv {{\sum\limits_{k = 0}^{2^{1 - 1_{- 1}}}\{ {{\Phi^{3D}\begin{pmatrix}{{X^{(1)}(k)};} \\{{Y^{(1)}(k)};} \\{\frac{B}{2^{1}},\sigma}\end{pmatrix}} + {\Phi^{3D}\begin{pmatrix}{{{{\hat{X}}^{({1 - 1})}(k)} - {X^{(1)}(k)}};} \\{{Y^{(1)}( {k + 2^{1 - 1}} )};} \\{\frac{B}{2^{1}},\sigma}\end{pmatrix}}} \}} + {2^{1}{\sum\limits_{k = 0}^{2^{1} - 1}{{\overset{\sim}{g}( \omega_{k}^{(1)} )}{{{\overset{\sim}{X}}^{(1)}( \omega_{k}^{(1)} )}}^{2}\mspace{14mu} {for}\mspace{14mu} \omega_{k}^{(1)}}}}}} = \frac{2\pi \; k}{2^{1}}},{0 \leq k < {2^{i}.}}} & (31)\end{matrix}$

Using the values determined in equations (29) and (30), the additionaloptimal values {circumflex over (X)}⁽¹⁾(τ+2¹⁻¹), 0≦τ<2¹⁻¹ are evaluatedas follows:

{circumflex over (X)} ⁽¹⁾(τ+2¹⁻¹)={circumflex over (X)}⁽¹⁻¹⁾(τ)−{circumflex over (X)} ⁽¹⁾(τ), 0≦τ≦2¹⁻¹  (32)

FIGS. 6A, 6B, 7A, and 7B show schematic images of phantoms. The imagesshown in FIGS. 6A and 7A were created using an ordered-subsetsexpectation maximization (OSEM) algorithm, as is known in the art. Theimages in FIGS. 6B and 7B were computed using equation (21), accordingto an embodiment of the present invention. For all images, phantoms wereconstructed by embedding six hollow or solid spheres, having diametersranging from 9.5 mm to 31.8 mm, in a cylinder. A first “hot spheres”phantom (FIGS. 7A, 7B) was constructed by filling the cylinder andhollow spheres with aqueous solutions of a radiochemical using ^(99m)Tc.The concentration of the solution in the spheres was twice theconcentration of the solution in the cylinder. A second “cold spheres”phantom (FIGS. 6A, 6B) was constructed in which solid non-radioactivespheres were used. For both phantoms, two-dimensional cross-sectionalimages were generated from one set of data taken by a camera having acollimator with parallel-sided channels.

Comparison of the images of FIGS. 6A, 7A with the images of FIGS. 6B, 7Billustrates a number of advantages of embodiments of the presentinvention:

-   -   The imaged radiochemical concentration in the cylinder, for both        types of phantom, is significantly more uniform in FIGS. 6B and        7B than in FIGS. 6A and 7A.    -   The ratio of radiochemical concentrations for the spheres and        for the cylinder is closer to the actual values in the images of        FIGS. 6B and 7B, compared to the image of FIGS. 6A and 7A.    -   The images of the spheres in FIGS. 6B and 7B are significantly        sharper than the corresponding images in FIGS. 6A and 7A.    -   The contrast between the imaged spheres and the cylinder in        FIGS. 6B and 7B is significantly greater than that between the        corresponding images in FIGS. 6A and 7A. Taken together with the        other differentiating factors described above, this leads to the        smallest sphere being very difficult to see in FIGS. 6A and 7A,        whereas it is clearly visible in FIGS. 6B and 7B.

It will be appreciated that the embodiments described above are cited byway of example, and that the present invention is not limited to whathas been particularly shown and described hereinabove. Rather, the scopeof the present invention includes both combinations and subcombinationsof the various features described hereinabove, as well as variations andmodifications thereof which would occur to persons skilled in the artupon reading the foregoing description and which are not disclosed inthe prior art.

1. A method for imaging, comprising: counting quanta of energy emittedinto a range of angles from particles of an energy emitter distributedover a volume in a body, thereby generating a set of counts; defining aprobability distribution expression that specifies a local concentrationof the particles of the energy emitter over the volume as a function ofthe set of counts, the function being defined in terms of respectivecoefficients of a plurality of different scales of the localconcentration, comprising at least a first and a second scale;specifying a dependence of the coefficients of the first scale on thecoefficients of the second scale; and computing the local concentrationover the volume by applying the probability distribution expression tothe set of counts subject to the specified dependence.
 2. The methodaccording to claim 1, wherein the dependence comprises a firstcardinality of the first scale.
 3. The method according to claim 2,wherein the dependence comprises a second cardinality of the secondscale.
 4. The method according to claim 1, wherein the first scale andthe second scale comprise sequential scales.
 5. The method according toclaim 4, wherein the dependence is given by:$\alpha_{s} = \frac{s^{2}}{( {s + n} )^{2}}$ wherein srepresents a cardinality of the first scale, 0.5≦n<1.5, and whereinα_(s) represents the dependence.
 6. The method according to claim 1,wherein an argument of the function comprises a term:(θ_(i) ^(s)−α_(s)θ_(i) ^(s+1)), wherein: i is an index of a location ofthe particles; s and s+1 are respective cardinalities of the first andsecond scales; θ_(i) ^(s),θ_(i) ^(s+1) represent respective coefficientsof the first and second scales; and α_(s) is a function of s.
 7. Themethod according to claim 1, wherein${f_{s}(z)} = \frac{{{z}_{2}}^{r}}{s^{r}\sigma}$ wherein: f_(s)(z)represents the function, having an argument z; s is a cardinality of thefirst scale; r is a positive number; and σ is a constant.
 8. The methodaccording to claim 7, whereinz=(θ_(i) ^(s)−α_(s)θ_(i) ^(s+1)), wherein: i is an index of a locationof the particles; s+1 is a cardinality of the second scale; θ_(i)^(s),θ_(i) ^(s+1) represent respective coefficients of the first andsecond scales; and α_(s) is a function of s.
 9. The method according toclaim 1, wherein the probability distribution expression comprises aterm−log [P(X)], wherein${{- {\log \lbrack {P(X)} \rbrack}} = {\sum\limits_{s}{\sum\limits_{i}{f_{s}( {\theta_{i}^{s} - {\alpha_{s}\theta_{i}^{s + 1}}} )}}}},$and wherein: X is a vector representing the local concentration;f_(s)(.) represents the function; i is an index of a location of theparticles; s and s+1 are respective cardinalities of the first andsecond scales; θ_(i) ^(s),θ_(i) ^(s+1) represent respective coefficientsof the first and second scales; and α_(s) is a function of s.
 10. Themethod according to claim 1, wherein the probability distributionexpression comprises a term−log [P(Y|X)], wherein${{- {\log \lbrack {P( {YX} )} \rbrack}} = {\frac{1}{2}{\sum\limits_{b = 1}^{b_{M}}\frac{( {y_{b} - ({HX})_{b}} )^{2}}{\Delta_{b}}}}},$and wherein: b is an index representing bins receiving the set ofcounts; b_(M) is a total number of the bins; y_(b) is a count receivedat a bin b; (HX)_(b) is an expected number of counts for the bin b Y isa vector representing the set of counts; X is a vector representing thelocal concentration; and Δ_(b) is a variance for the bin b.
 11. Themethod according to claim 10, whereinΔ_(b)=max(y _(b) ,B), wherein B is a constant.
 12. The method accordingto claim 1, wherein the probability distribution expression comprises aterm having a distribution chosen from a Poisson distribution, aLaplacian distribution, and a non-white Gaussian distribution.
 13. Themethod according to claim 1, wherein the coefficients comprise waveletcoefficients of a Haar wavelet.
 14. The method according to claim 1,wherein computing the local concentration comprises performing aniteration on the probability distribution expression to find a maximumvalue of a probability distribution comprised in the expression.
 15. Themethod according to claim 1, wherein the local concentration comprises aplurality of time-dependent-local-concentrations comprising a firsttime-dependent-local-concentration measured at a region of the volumeduring a first time slot and a secondtime-dependent-local-concentration, different from the firsttime-dependent-local-concentration, measured at the region during asecond time slot, and wherein the probability distribution expressioncomprises a time-dependent regularization function relating the firsttime-dependent-local-concentration and the secondtime-dependent-local-concentration.
 16. The method according to claim15, wherein the time-dependent regularization function comprises aFourier transform of a sequence of time-dependent-local-concentrationsat the region, the sequence comprising the firsttime-dependent-local-concentration and the secondtime-dependent-local-concentration.
 17. The method according to claim15, and comprising synchronizing at least one of the first time slot andthe second time slot in response to a time-dependent signal receivedfrom the body.
 18. The method according to claim 15, wherein computingthe time-dependent-local-concentrations comprises iteratively computingpartial sums of the time-dependent-local-concentrations.
 19. The methodaccording to claim 1, wherein the plurality of different scalescomprises a number of scales chosen from between three and seven scales.20. The method according to claim 19, wherein the number of scalescomprises five scales.
 21. The method according to claim 1, wherein theprobability distribution expression comprises a conditional expectation.22. Imaging apparatus, comprising: a camera which is arranged to countquanta of energy emitted into a range of angles from particles of anenergy emitter distributed over a volume in a body, thereby generating aset of counts; and a processing unit which is configured to: define aprobability distribution expression that specifies a local concentrationof the particles of the energy emitter over the volume as a function ofthe set of counts, the function being defined in terms of respectivecoefficients of a plurality of different scales of the localconcentration, comprising at least a first and a second scale, specify adependence of the coefficients of the first scale on the coefficients ofthe second scale, and compute the local concentration over the volume byapplying the probability distribution to the set of counts subject tothe specified dependence.
 23. The apparatus according to claim 22,wherein the local concentration comprises a plurality oftime-dependent-local-concentrations comprising a firsttime-dependent-local-concentration measured at a region of the volumeduring a first time slot and a secondtime-dependent-local-concentration, different from the firsttime-dependent-local-concentration, measured at the region during asecond time slot, the apparatus comprising: a detector which detects atime-dependent signal from the body and which conveys the time-dependentsignal to the processing unit so as to synchronize at least one of thefirst time slot and the second time slot.
 24. A computer softwareproduct for imaging, comprising a computer-readable medium havingcomputer program instructions recorded therein, which instructions, whenread by a computer, cause the computer to: receive counts of quanta ofenergy emitted into a range of angles from particles of an energyemitter distributed over a volume in a body, thereby generating a set ofcounts, define a probability distribution expression that specifies alocal concentration of the particles of the energy emitter over thevolume as a function of the set of counts, the function being defined interms of respective coefficients of a plurality of different scales ofthe local concentration, comprising at least a first and a second scale,specify a dependence of the coefficients of the first scale on thecoefficients of the second scale, and compute the local concentrationover the volume by applying the probability distribution expression tothe set of counts subject to the specified dependence.